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Abstract 

We present a detailed derivation and thin interface analysis of a phase-field model that can accu- 
rately simulate microstructural pattern formation for low-speed directional solidification of a dilute 
binary alloy. This advance with respect to previous phase-field models is achieved by the addition 
of a phenomenological "antitrapping" solute current in the mass conservation relation [A. Karma, 
Phys. Rev. Lett 87, 115701 (2001)]. This antitrapping current counterbalances the physical, albeit 
artificially large, solute trapping effect generated when a mesoscopic interface thickness is used to 
simulate the interface evolution on experimental length and time scales. Furthermore, it provides 
additional freedom in the model to suppress other spurious eff'ects that scale with this thickness 
when the diflFusivity is unequal in solid and liquid [ R. F. Almgren, SIAM J. Appl. Math 59, 2086 
(1999)], which include surface diffusion and a curvature correction to the Stefan condition. This 
freedom can also be exploited to make the kinetic undercooling of the interface arbitrarily small 
even for mesoscopic values of both the interface thickness and the phase-field relaxation time, as 
for the solidification of pure melts [A. Karma and W.-J. Rappel, Phys. Rev. E 53, R3017 (1996)]. 
The performance of the model is demonstrated by calculating accurately for the first time within 
a phase-field approach the MuUins-Sekerka stability spectrum of a planar interface and nonlinear 
cellular shapes for realistic alloy parameters and growth conditions. 
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I. INTRODUCTION AND SUMMARY 



In recent years, the phase-field method has become a standard tool to simulate microsc- 
tructure evolution in materials [l], a subject of both fundamental and applied interest 
and more generally to tackle free boundary problems. Its chief advantage is to avoid front 
tracking by making phase boundaries spatially diffuse with the help of order parameters, 
termed phase fields, which vary smoothly between bulk phases. 

Simulating the evolution of complex morphologies in two and three dimensions is in prin- 
ciple straightforward with this method. Making quantitative predictions on experimentally 
relevant length and time scales, however, has been a major challenge. This challenge stems 
from the fact that phase-field simulations are simply not feasible if parameters of the model 
are chosen to match those of a real physical system. With such a choice, both the width W of 
the diffuse interface and the characteristic dissipation time scale r in the phase-field kinetics 
are microscopic: is a few Angstroms and r is roughly the ratio of W and the thermal 
velocity of atoms in the liquid jj, ^ 0]. In contrast, diffusive transport of solute in bulk 
phases occurs on macroscopic length and time scales that are several orders of magnitude 
larger than W and r, respectively. Therefore, resolving both microscopic and macroscopic 
length/time scales in phase-field simulations for typical experimental solidification rates of 
fim/sec to mm/sec is impractical, even with efficient algorithms. 

In view of this, the only possible choice is to carry out simulations with W and r orders 
of magnitude larger than in a real material. The question becomes then whether the phase- 
field model is still quantitatively meaningful with such a choice. The rest of this section 
explores the answer to this question in the context of previous works and serves both as a 
summary and a guide for the following sections of this paper. To conclude this section, we 
summarize the mains results needed to carry out quantitative simulations of the directional 
solidification of a dilute binary alloy. 



A. Capillarity 

In the phase-field model of a pure substance (of say A molecules) , the excess free-energy of 
the solid-liquid interface, 7, is determined by the combination of the bulk free-energy density 
at the melting point, T^), which is a double- well function with minima corresponding to 
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solid and liquid and the gradient square term, (t| V0p. Minimization of the total free-energy, 
which is the spatial integral of the sum of these two terms, yields the standard result that 
7 ~ WH, where H is the barrier height of the double-well potential, and W ~ {a/HY^'^ 
is the width of the tanh-profile in the diffuse interface. This results implies that there 
always exist a pair of values of a and H for any pair of values of W and 7. Thus, the 
experimental magnitude of 7 in the classic Gibbs-Thomson condition can be reproduced even 
if a computationally tractable "mesoscopic" interface thickness (i.e. on a scale comparable 
to the microstructure) is used in the phase-field model. Optimally, this thickness should be 
chosen just small enough to resolve accurately the interface curvature. 

A phase-field model for a dilute alloy can generally be constructed by adding to the 
free-energy density the contribution of solute B molecules. The simplest way to construct 
this free-energy is to interpolate between the known free-energy densities in solid and liquid 
with a single function of 0, as in the original model of Wheeler et al. j^] (see also Ref. ^). 
From a computational standpoint, however, this approach places a stringent constraint on 
the interface thickness. The reason is that there is generally an extra contribution to 7 due 
to solute addition that depends on interface thickness, solute concentration at the interface, 
and temperature. In section III. A, we show how this extra contribution can generally be 
made to vanish by using two different functions of 0, which interpolate separately between 
solid and liquid the enthalpic (internal energy) and entropic part of the free-energy density. 
The condition that this contribution vanishes takes the form of an algebraic relation between 
these two interpolation functions. If this relation is satisfied, the model introduced previously 
in Ref. P] is recovered. The equilibrium phase-field profile decouples from the equilibrium 
solute concentration profile and 7 ~ WH, as for a pure substance. This removes the 
constraint on the interface thickness associated with solute addition without the need to 
introduce separate concentration fields in each phase as in Refs. 
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B. Interface-thickness-dependent nonequilibrium effects 

The main conclusion from the preceding paragraphs is that the phase-field method pro- 
vides sufficient freedom to choose W arbitrarily large to model capillarity. However, mi- 
crostructural pattern formation is also generally controlled by nonequilibrium effects at the 
interface. For a microscopic W and low solidification velocities, these effects are negligibly 
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small. The interface relaxes rapidly to a local thermodynamic equilibrium and its non- 
linear evolution is driven by slowly evolving gradients of thermodynamic quantities in bulk 
phases. For a mesoscopic thickness, however, these nonequilibrium effects become artificially 
magnified, thereby competing with, or even superseding, capillary effects, and dramatically 
altering the large scale pattern evolution. Therefore, the central challenge of quantitative 
phase-field modeling of solidification at low velocity, onto which we focus in the present 
work, consists of formulating the model, and knowing how to choose its parameters, in order 
to avoid unphysically large non-equilibrium effects at the interface. This is in contrast to 
rapid solidification where nonequilibrium effects play a dominant role. In this case, the chal- 
lenge consists of describing the correct magnitude of these effects with mesoscale phase-field 
parameters, which requires a different approach (see Ref. 12|). 

For pure materials. Karma and Rappel ^ have developed a thin interface analysis, which 
only assumes that W is small compared to the scale of the microstructure. This analysis 
shows that the standard free-boundary problem of solidification — a classic Stefan condition 
together with a velocity-dependent form of the Gibbs-Thomson relation that incorporates 
interface kinetics — is recovered even for a mesoscopic W. Heat diffusion in the mesoscale 
interface region only modifies the expression for the interface kinetic coefficient, /i^. This 
"renormalization" of /i^ has the crucial property that r needs not be microscopic to make 
this coefficient arbitrarily large (arbitrarily fast kinetics), and hence to simulate the limit of 
local equilibrium at the interface dominated by capillarity. 

This advance bridges the gap between the atomistic scale of interfacial phenomena and the 
mesoscale of the microstructure. In addition, efficient multi-scale simulation algorithms have 
been developed to bridge the remaining gap between the microstructure and the transport 
scales ^ 3| • The combination of these two advances has lead to the first direct quantitative 
comparison between fully three-dimensional phase-field simulations of dendritic growth in 
pure melts at low undercooling and experiments [lfi| . 

Achieving the same success for alloys has turned out to be considerably more difficult. 
A major source of difficulty is that solute diffusion is generally much slower in solid than 
liquid. When diffusion is asymmetrical, the use of a mesoscopic W artificially magnifies 
several nonequilibrium effects at the interface that are absent when diffusion is symmetrical. 
Co nsequen tly, phase-field models in which one or several of these effects are present 0, Isl 
are not suitable for quantitative simulations at low velocity. 



Consequen 
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These nonequilibrium effects were first cliaracterized in detail by Almgren using a tliin 
interface analysis of a phase-field model of the solidification of pure melts with asymmetric 
diffusion. Directly analogous effects are present in alloy solidification which include (i) 
solute diffusion along the arclength of the interface (surface diffusion), (ii) a modification 
of mass conservation associated with the local increase of arclength of a moving curved 
interface (interface stretching), and (iii) a discontinuity of the chemical potential of the 
dilute impurity across the interface. 

These nonequilibrium effects originate physically from solute transport in the mesoscale 
interface region that is governed by the standard continuity equation for a dilute alloy 

dc Dvq 



-V ■ (g(0)cV/i) , (1) 



dt RT„^ 

where R is the gas constant, Vq is the molar volume of solute molecules, is the melting 
temperature, /i is the chemical potential, and the product Dq{(j)) governs how the solute 
diffusivity varies through the diffuse interface, from zero in the solid (for a one-sided model) 
o a constant value D in the liquid. The best known of these effects is solute trapping 



3, 



2l| that is associated with the chemical potential jump at the interface. The problem 



is that the magnitude of all these effects scales with the interface thickness. Since W in 
phase-field computations is orders of magnitude larger than in reality, solute trapping will 
become important at growth speeds where it is completely negligible in a real material. 
Surface diffusion and interface stretching, in turn, modify the mass conservation condition 

ci{l-k)Vr, = -D— + .-- (2) 

where q is the concentration on the liquid side of the interface, k is the partition coefficient, 
Vn is the normal interface velocity, and " ■ ■ ■ " is the sum of a correction ~ q(1 — k)WVn}C, 
corresponding to interface stretching, where JC is the local interface curvature, a correction 
~ WDd'^ci/ds'^ , corresponding to surface diffusion along the arclength s of the interface, and 
a correction ~ kci{l—k)WV^ / D proportional to the chemical potential jump at the interface. 
All three corrections, which are proportional to the interface thickness, are negligible in a 
real material at low velocity. For this reason, they have not been traditionally considered in 
sharp- interface models (reviewed in Sec. II). For a mesoscopic interface thickness, however, 
the magnitude of these corrections becomes comparable to the magnitude of the normal 
gradient of solute, thereby modifying Vn and the pattern evolution. Thus, the phase-field 
model must be formulated to make all three effects vanish. 
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C. Limitation of variational models 



The model discussed in Sec. III. A follows the general approach of nonequilibrium ther- 
modynamics where the evolution equations for cj) and c are derived variationally from a 
Lyapounov functional T that represents the total free-energy of the system. The resulting 
"gradient dynamics" guarantees that decreases monotonously in time in an isolated sys- 
tem. In addition to the double- well potential /(</>), this variational model contains three 
basic interpolation functions: the two functions that interpolate between solid and liquid 
the enthalpic and entropic part of the free-energy density (section III. A), which we denote 
here by ^(0) and ^(0), respectively, and the diffusivity function g(0) in Eq. (jT)) that varies 
from zero in solid to unity in liquid. 

These functions should in principle be chosen to cancel all spurious interface-thickness 
dependent effects. As already discussed in Sec. LA, a quantitative description of capillarity 
can be obtained by requiring that the solute contribution to 7 vanish. This condition is 
only satisfied if the two functions (7(0) and ^(0) are related, and the latter determines the 
equilibrium solute concentration profile 

Co(0) = ^ + ^(0)^, (3) 

where ^(0) varies from +1 in the solid where Cq = Cg to —1 in the liquid where Cq = Q. 

We are left with only two functions, ^(0) and g(0), to satisfy the three aforementioned 
conditions that surface diffusion, interface stretching, and the chemical potential jump at the 
interface, should vanish. The thin-interface analysis of section IV applied to this variational 
model shows that these three conditions are given, respectively, by 

fO r+00 

dr q{(j){r))cQ{(j){r)) = / rfr [q - g(0(r))co(0(r))] , (4) 

00 JO 

r+00 

rfr [co(0(r)) - c^] = / rfr [q - co(0(r))] , (5) 



00 JO 



00 g(0(r))co(0(r)) Jo 



-00 



co(0(r)) - Cs 



q{(l){r))coi(/){r)) 



(6) 



where k = Cg/ci is the partition coefficient, r is the coordinate normal to the solid-liquid 
interface that varies from —00 in solid to +00 in liquid far from the interface, and Cq is given 
by Eq. (jH)) that can be assumed to remain valid for a slowly moving interface. 

A simple physical interpretation of these conditions is obtained by analogy with Gibbs' 
treatment of interfacial phenomena where "excess quantities" are attributed to a mathemat- 




FIG. 1: Illustration of the definition of surface excess. The excess of solute is the integral along r 
of the actual solute profile (thick solid line) minus its step profile idealization (thick dashed line) 
with the Gibbs dividing surface at r = 0. This excess is negative in the depicted example. The thin 
solid line depicts the phase- field profile, (l){r) = — tanh(r/\/2VF). The standard mass conservation 
condition [Eq. is recovered if all three excess quantities defined by the difference between the 
left-hand-side and the right-hand-side of Eqs. Q-®, vanish. 

ical surface with zero volume dividing two phases, which corresponds here to r = 0. In this 
analogy, Eqs. are conditions that excess quantities of the interface vanish. For exam- 

ple, as illustrated in Fig. [H the excess of solute is the integral through the diffuse interface 
of the difference between the actual smoothly varying solute profile cq and the imaginary 
step function profile equal to Cs for r < and q for r > 0. The condition that this excess 
vanishes is identical to Eq. It implies that mass conservation is left unchanged if there 
is no excess of solute to redistribute along the arclength of the interface. Similarly, surface 
diffusion vanishes [Eq. if there is no excess of the transport coefficient q{4>)c multiplying 
the chemical potential gradient in Eq. (P). Finally, the jump of chemical potential vanishes 
if there is no excess of chemical potential gradient [Eq. ©]. This condition is simple to 
derive for a flat interface by rewriting Eq. in a local frame moving at velocity V (i.e. 
d/dt —>■ —Vd/dr and V d/dr). After integrating both sides of Eq. (0) once with re- 
spect to r, one obtains the expression for the chemical potential gradient through the diffuse 
interface dfi/dr ^ —V{co — Cs)RTm/ (Dvoqco) , and hence Eq. (0). 

A major pitfall of this variational model is that all three excess quantities cannot be 
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made to vanish simultaneously with only two functions g{(f)) and q{(f)). For example, with 
the standard quartic form of the double-well, which is an even function of (p, the equilibrium 
profile is an odd function of r. Therefore, Eq. (0) can be made to vanish by choosing ^((/)) 
to be an odd function of (p. It is then possible to choose q{(f)) to satisfy Eq. (@)). However, 
this leaves no freedom to make the jump of chemical potential vanish. More generally, it is 
possible to make two of the three excess quantities vanish for different choices of g{(f)) and 

but not the three of them simultaneously. 

1 I 

Elder et al. [22] proposed to make the discontinuity of chemical potential vanish by 
an appropriate choice of interface position (Gibbs dividing surface) which makes the cor- 
responding excess quantity vanish. These authors, however, did not take into account the 
other two excess quantities found by Almgren for asymmetric diffusion [18] . These quanti- 
ties appear at higher orders in the asymptotic expansion used by Elder et al. which, for the 
solidification of pure melts with symmetrical diffusion, yields the same results as the thin 
interface analysis of Karma and Rappel For asymmetrical diffusion, all three excess 

quantities can generally not be made to vanish by a redefinition of the interface position. 

It might be possible to make all three excess quantities vanish for non-trivial oscillatory 
forms of the functions ^(0) and Such forms, if they exist, would require a high res- 

olution of the interfacial layer that is not computationally desired. Also, other variational 
models than the one discussed here are in principle possible. McFadden et al. have formu- 



lated a variationa 
conductivities l2 



phase-field model of the solidification of pure melts with unequal thermal 
. This model provides additional freedom to cancel the discontinuity of 
temperature at the interface, which they interpret as "heat trapping" by analogy with solute 
trapping that is associated with the discontinuity of chemical potential in the case of alloys. 
However, as Elder et al., these authors did not consider the additional constraints associated 
with surface diffusion and interface stretching for a non-planar interface. While we cannot 
rule out that it may be possible to formulate variational models that remove all constraints 
on the interface thickness, achieving this goal appears extremely difficult. 



D. Non-variational models and antitrapping 

A way out of this impasse is to drop the requirement that the equations of the phase- 
field model be strictly variational. This provides additional freedom to cancel all spurious 
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corrections produced by a mesoscale interface thickness. As shown recently in Ref. jioj ]. 
a successful approach consists of adding a phenomenological "antitrapping current" in the 
continuity relation [Eq. This current produces a net solute flux from solid to liquid 

proportional to the interface velocity that counteracts solute trapping and restores chemical 
equilibrium at the interface. By adjusting the magnitude of this current, which modifies Eq. 
it is therefore possible to satisfy simultaneously Eqs. (jH)-®- 
Furthermore, the same function g{(f)) must appear in the evolution equations for (j) and the 
continuity relation [Eq. (Q)] in the variational model. The additional freedom to replace ^(0) 
by another function h{(f)) in the modified continuity relation with the antitrapping current 
turns out to be critically important to obtain the same renormalization of the interface 
kinetic coefficient as in the analysis of Karma and Rappel for pure melts jisl . 



E. Summary of phase-field equations and thin-interface limit 

We summarize here the equations of the non-variational phase-field model for the di- 
rectional solidification of a dilute binary alloy that are needed to carry out quantitative 
computations. The lengthy details of the derivation of the model and of the asymptotic 
analysis are exposed in sections III and IV below. The model uses the standard low velocity 
frozen temperature approximation, T = Tq + G{z — Vpt), where Vp is the pulling speed and 
G is the temperature gradient. The basic equations of the model are 

. W'V^, + , - ,3 _ (e- - 1 - ^) . (7) 

^ = v(£>9WcVm-j.,), (8) 



where 



is a dimensionless measure of the deviation of the chemical potential from its equilibrium 
value at a reference temperature Tq with corresponding liquidus concentration c°, m < 
is the liquidus slope, 

= -aW{l - k)cy^^, (10) 
^ ^ ' dt |V0| 

is the antitrapping current, 

^(T) = ro(l + (11) 



mcr 
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is a temperature-dependent phase-field relaxation time, and A is a dimensionless coupling 
constant. For the choices h{(j)) = 0, g(0) = (1 - 0)/[l + A; - (1 - k)(f)], g{(j)) = (15/8)(0 - 
20^/3+0^/5), and a = l/{2^/2), this model reduces in its thin-interface limit to the standard 
one-sided model of alloy solidification. The chemical capillary length and the interface 
kinetic coefficient (3 (defined in section II) are related to the phase-field parameters by 

do = aiW/X, (12) 



02 



ToD 



(13) 



where A = 15A/8 and ai = 5-\/2/8 and 02 = 0.6267 are the same numerical constants as in 
Ref. we note that A has been defined for convenience in the present paper to avoid 

carrying a numerical factor of 15/8 in the thin- interface analysis of the equations. 

A previous version of this model for isothermal alloy solidification was presented in Ref. 
[igj l together with benchmark computations for dendrite growth. The present extension 
to non-isothermal growth conditions introduces a temperature-dependent relaxation time 
r(T). As discussed in more details in section IV. C, this new feature makes it possible to 
achieve vanishing interface kinetics (i.e. local equilibrium at the interface) for the entire 
range of interface temperature that occurs during directional solidification. For simplicity, 
we have written down the equations of the model for isotropic surface tension and interface 
kinetics. The extension to anisotropic growth is discussed in section IV.E. Also, both for 
simulating and analyzing the above equations, it is convenient to rewrite them in terms of a 
new variable U = (e" — 1)/(1 — k). This avoids numerical computations of exponential and 
logarithm functions. In addition, it transforms the equations in a form closely related to the 
phase-field model for the solidification of a pure substance where U is the direct analog of 
the temperature field. Details of this change of variable are given in section III.B.c. 

Simulations of microstructural pattern formation using this model are presented in section 
V, which also contains the final form of the anisotropic phase equations ()132|) and ()133p that 
are solved numerically. We report the first ever quantitative phase-fie^ 
classic MuUins-Sekerka linear stability spectrum of a planar interface 



d computation of the 
and nonlinear cell 



shapes for realistic experimental parameters of low velocity directional solidification. 
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II. SHARP-INTERFACE MODELS 



We consider the solidification of a dilute binary alloy made of substances A and B, with 
an idealized phase diagram that consists of straight liquidus and solidus lines of slopes m 
and m/k, respectively, where k is the partition coefficient. The interface is supposed to be 
in local equilibrium, that is, 

Cs = kci, (14) 

where Cg and q are the concentrations (in molar fractions) of impurities B at the solid and 
liquid side of the interface, respectively. 

The interface temperature satisfies the generalized Gibbs- Thomson relation, 

T = Tm-\m\ci-TIC-Vn/fJ,k, (15) 
where Tm is the melting temperature of pure A, 

r = (16) 

the Gibbs- Thomson constant, 7, the surface tension, L, the latent heat of fusion per unit 
volume, /C, the interface curvature, Vn its normal velocity, and fik the linear kinetic coeffi- 
cient. Here, the surface tension and the kinetic coefficient are taken isotropic for simplicity; 
anisotropic interface properties will be considered below. 

Heat is supposed to diffuse much faster than impurities, so that the temperature field 
can be taken as fixed by external conditions, in spite of the rejection of latent heat during 
solidification. Then, Eq. (P3j) yields a boundary condition for the solute concentration at 
the interface. 

Of particular interest is the one-sided model of solidification that assumes zero diffusivity 
in the solid. This is often a good approximation for alloy solidification, in which the solute 
diffusivity in the solid may be several orders of magnitude lower than in the liquid. 



A. Isothermal solidification 

For isothermal solidification at a fixed temperature Tq < Tm, the concentration obeys the 
set of sharp-interface equations 

dtc = DV^c, (17) 
12 



q(1 - k)Vn = -DdnC\ + , (18) 

q/cO = 1 - (1 - k)dolC~ (1 - k)(3Vn (19) 

where D is the solute diffusivity in the hquid, Ki, the normal velocity of the interface, dnc\^, 
the derivative of the concentration field normal to the interface, taken on the liquid side of 
the interface, 

c° = (T„-To)/|m|, (20) 



the equilibrium concentration of the liquid at Tf 



0, 



4 = ^, (21) 

the chemical capillary length, where ATq = |m|(l — k)c^ is the freezing range, and f3 = 
l/[/ifcATo]. Equation (fTH|) . the Stefan condition, expresses mass conservation; Eq. (fTTIjl can 
be directly obtained from Eq. (fTH|) . 

B. Directional solidification 

For directional solidification, we use the frozen temperature approximation, in which the 
temperature field for solidification with speed in a temperature gradient of magnitude G 
directed along the z axis is taken as 

T{z)=To + G{z-Vpt). (22) 

Now To is given by inverting Eq. (|^. and cj* = Coo/k, where Coo = c{z = +oo) is the global 
sample composition. Thus, is the solute concentration on the liquid side of a steady-state 
planar interface. Then, Eq. ()19p is replaced by 

q/c° = 1 - (1 - k)do /C - (1 - k)f3Vn -{l-k){z- Vpt)/lT (23) 

where 

, _ (24) 



a 

is the thermal length. 
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C. Formulation in terms of dimensionless supersaturation 



In order to later compare with the sharp-interface hmit of the phase-field models treated 
here, we rewrite Eqs. (jl7ll8l23j) in terms of the local supersaturation with respect to the 
point (c°,To), measured in units of the equilibrium concentration gap at that point, 



c-q 



,0 



We obtain 



dtU = DV'^U (liquid), (26) 
[l + {l-k)U]Vn = -DdnU\+ (interface), (27) 
U = -doJC- pVn-iz-Vpt)/lT (interface). (28) 

Note that, for k = 1, we recover the constant miscibility gap model. Furthermore, if we 
reinterpret f/ as a dimensionless temperature and drop the directional solidification term 
{z — Vpt)/lT, we obtain a one-sided version of the pure substance model. 

III. PHASE-FIELD MODELS 

In this section we first derive a generic variational model (Sec. IIII A|) . and we then modify 
it in view of canceling spurious effects (Sec. IIIIB|) . 

A. Variational formulations 

In a phase-field model, a continuous scalar field (p is introduced to distinguish between 
solid (0 = +1) and liquid (0 = —1). The two-phase system is usually described by a 
phenomenological free energy functional. 



F[0,c,T] 

where 



dV 



-|V0P + /(0,T„) + /^b(0,c,T) 



(29) 



/(0,T^) = iJ(-0V2 + 0V4) (30) 



is the standard form of a double-well potential providing the stability of the two phases 
= ±1 with a barrier height H, /yiB(0, c, T) changes their relative stability as a function 

14 



of the position in a T-c phase diagram, and the term in a provides a penalty for phase 
gradients which ensures a finite interface thickness. H has dimensions of energy per unit 
volume, and a of energy per unit length. 

In a variational formulation, the equations of motion for all fields (here the concentration 
and phase fields) can be derived from that functional: 

|^V.(M(0,c)vf), (32) 

where K^{T) is a kinetic constant that can generally be temperature-dependent. The second 
equation is a statement of mass conservation, since it can be rewritten as 

dc -> -> 

— + V ■ Jc = 0, (33) 

— » — * 

where Jc = —MV^i is the solute current density, /i = 5F/5c is the chemical potential, and 
M(0, c) is the mobility of solute atoms or molecules, which we choose to be 

M(</.,c) = ^Dg(0)c (34) 

in order to later obtain Pick's law of diffusion in the liquid. Here, Vq is the molar volume 
of A, R, the gas constant, and q{4>), a dimensionless function that interpolates between in 
the solid and 1 in the liquid, and hence dictates how the solute diffusivity varies through the 
diffuse interface. Note that we have not included an equation of motion for the temperature 
field, since we consider it fixed by external constraints. Of course, the formalism could be 
extended to include an appropriate equation for heat transfer 2^. 

An important step is the construction of the function fj^B that interpolates between the 
free energy densities of the bulk phases (solid and liquid). While these bulk free energies 
should reduce to the curves that can be obtained from thermodynamic databases, the de- 
pendence of fyiB on (f) influences only the interfacial region, and this freedom can be used 
to construct a particularly simple phase-field model. This will be illustrated here for the 
case of a dilute binary alloy. First, we consider the bulk free energies and make sure that 
they reproduce the equilibrium properties of the sharp-interface model of Sec. IHl Then, we 
interpolate between them. 
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For a dilute alloy, the free energies of solid and liquid /jy(c, T) can be written as the sum 
of the free energy of pure A, f^{T), and contributions due to solute addition: 

U{c,T) = f^{T) + {c\nc-c)+e,c z/ = 1, s. (35) 

The second term on the right hand side is the dilute form of the mixing entropy, and the 
term EiyC is the change of the internal energy density. We expand this expression to first 
order in T — Tm to recover the straight liquidus and solidus lines of Sec. ITTl 

T) = f^{T^) - s,{T - T^) + ^(clnc - c) + e,c, (36) 

where Si, = —df^/dT are the entropy densities of solid and liquid at T^. By using 
instead of T in the mixing entropy, we have neglected terms of order (T — Tm)c, which are 
second-order for dilute alloys. 

The phase diagram is determined by the standard common tangent construction, which 
is equivalent to requiring that the chemical potential and the grand potential u (i.e., the 
thermodynamic potential for a varying number of solute particles) be equal in the solid and 
liquid. The corresponding equilibrium concentrations Cs{T) and q(T) are the solutions of 



dfs{c,T) 



dc 



dfiic,T) 



dc 



= f^E{T), (37) 

C=Cl 

fs{Cs,T) - fiECs = fl{ci,T) - fiECl = UJe(T). (38) 



The first equality yields the partition relation Eq. (fT^ . = kci, with a partition coeffi- 
cient 

.^exp(-||£). (39) 



where we have defined Ae = Eg — Si. Combining this result with Eq. (jH8|l yields 

(40) 

where we have used that the latent heat per unit volume is L = Tm{si — Ss). From Eq. (HOJ, 
we identify the liquidus slope to be 

,n^-mifJi>, (41) 

VqL 

the van't Hoff relation for dilute binary alloys. 
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In the standard phase-field approach, the two bulk free energies are interpolated with the 
help of a single function of the phase field 0. Here, it is advantageous to use two different 
interpolation functions for the entropy and the internal energy terms, 

RT 

UB{<f>, c, T) = /^(T^) - (T - Trn)s{<f)) + ^{c\nc-c)+ e{<f))c, (42) 

with 

5(0) = 6 + g{(j)) Ae/2, (43) 
^(0) = ^-^(0)^, (44) 

where e = {es + ei)/2, and we have again used L = Tm{si — Ss) in s(0). ^(±1) = ^?(±1) = ±1, 
and we further require ^'(±1) = g'{±l) = for = ±1 to remain bulk equilibrium solutions 
for any value of c and T. 

This completes the model specification, except for the interpolation functions g{(j)) and 
g{4>). In order to choose them appropriately, it is important to consider the equilibrium 
properties of the model, which follow from the conditions 

-J- = l^E, (45) 

where /x^; is the spatially uniform equilibrium value of the chemical potential. These two 
equations uniquely determine the spatially varying stationary profiles of c and in the diffuse 
interface region, Co(x) and (f)o{x). Since the phase field interpolates between the two bulk free 
energies, the limiting values of the concentrations and the equilibrium chemical potential are 
the ones determined by the common tangent construction above. From Eq. (03)), we have 

RTm - -/I \ A^ / An\ 

Inco + £ + g[(t)o)— = fJ'E, (47) 

Vo I 

from which we obtain the expression for the equilibrium concentration profile using the 
solution of Eq. ^ and Eq. (jSHI), 

Co(a;) = Q exp [1 + M0o(a:))]) = c,k'^^^-^^'^^^-))\l\ (48) 

From the equilibrium condition for 0, Eq. we obtain 
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With the help of Eqs. (j^ . (jlUI) and (jHJ, the right-hand side can be rewritten as 



(l-A;)^'(0o)+lnA;^^'(0o) 



(50) 



For a generic choice of the functions g and g, and in particular for the "standard" choice 
g = g, no analytic solution for is known. Furthermore, the equilibrium solution and its 
properties, in particular its surface tension, depend on the various coefficients that appear in 
the right-hand side. This can be avoided if the right-hand side vanishes {d^fAB{<Po, co, T) = 
0). With the help of Eq. ()48|). we obtain the condition on the interpolation functions, 

(l_,)iM + ln*iMexp(!^|l+SW|)=0, (51) 

It can be used to eliminate one of them in terms of the other. Taking into account the 
requirement ^(±1) = '^(±1) = ±1, we find 

l + fc_2exp(i^ [1 + ^(0)]) 1 + A; - 2A;[i+5Wl/2 ^ ^ 

^('^^ = r^k = r^k ' 

Using the latter relation, the equilibrium concentration profile can also be rewritten as 

co(0) = Q ^ = + gW— — • (54) 

The physical meaning of the two interpolation functions is hence completely transparent: g 
interpolates the internal energy [Eq. ()43p]. and as a consequence the chemical potentials [Eqs. 
(jTTj) and PH|) ]. whereas g interpolates the entropy density [Eq. (jH))] and, as a consequence 
of Eq. (jnH), the concentration [Eq. (jSlj)]- 

If Eq. (15111 is satisfied, the right-hand side of Eq. ()5()j) vanishes, and the solution for the 
equilibrium profile of is the usual hyperbolic tangent, (poix) = — taioh{x / ^/2W) , where 
W = (a/Hy^"^ measures the width of the diffuse interface. Furthermore, the surface tension 
is defined as the excess of the grand potential = / — fic, integrated through the interface, 
that is, 7 = J^^Lj{x) —uje- Because condition (jKTjl is equivalent to require d(i,fAB{4>o, cq, T) = 
0, under this condition fAB{<Po,co,T) is independent of x and equals its bulk phase values 
fy{cu,T). Since the latter enter the expression for the equilibrium grand potential uje as 
given by Eq. (jHHj) . the contribution of fAs to uj{x) — uje is zero. Thus, only the two other 
interface terms in Eq. contribute. Taking into account that both contribute the same 
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amount (equipartition relation), we have u^x) — ue = H\l — 0o(a;)^]^/2, and hence the 
surface tension is 

7 = IWH (55) 

with / = 2-\/2/3. Like in the sharp-interface model of Sec. m 7 is independent of solute 
concentration and temperature. Let us stress again that this property is only achieved if 
condition (jKTjl is satisfied. Otherwise, Eq. (jKKjl is replaced by a more complicated expression 
which contains the impurity concentration, and which needs in general to be calculated nu- 
merically. A drawback of this more complicated expression is that the dependence of 7 on 
concentration along the interface cannot be chosen independently of the value of W . This 
feature leads to an unphysically large variation of 7 with concentration for computationally 
tractable mesoscopic values of W . Eq. (j55|) yields a concentration-independent expression 
for 7 that is free of this limitation. Moreover, the fact that the equilibrium profile remains 
a hyperbolic tangent for arbitrary values of the concentration makes the relationship be- 
tween phase-field and sharp-interface parameters obtained from the thin-interface analysis 
independent of the value of the local concentration. This, in turn, avoids spurious kinetic 
corrections that are present otherwise. 

Once we have found a convenient relation between ^(0) and ^(0), we come back to the 
complete dynamical model. The relations we have found in equilibrium can now be used to 
obtain two particularly simple forms of the phase-field equation out of equilibrium. For the 
first, we remark that Eq. ()51|) implies that g'{(f)o)ci{l — k) = ~g'{(j)Q)lnkcQ, and therefore 
the function g can be eliminated in favor of the phase-dependent equilibrium concentration 
Co(0, T) and the function g. Dividing Eq. (jHT|) by H, we obtain 

r— = iyV + 0-0 + — — In/c^ 

at 2voHm 

with r = l/[Kfi,{T)H]; the driving force is the local supersaturation. The temperature 
dependence of r will be addressed later in section IV. C. 

The second possibility is to rewrite the phase-field equation in terms of the dimensionless 
variable 

" = - = (^/^°) - ¥^ ^'^^"^^ + = ( c°[l + A:-(l-fc)^(0)] ) ' ^'^^ 

which measures the departure of the chemical potential from its value ^e{To) for a flat 
interface at the equilibrium liquidus temperature Tq (and liquid concentration = q(To)). 
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c- co(0,T) 
ci{T) 



(56) 



Then, it is preferable to eliminate g{(f)) in favor of g{(f)). The result is the form used in 
Ref. Iia, 

- + ^ - - t4i^'(*) - 1 - ^) • 

where we have defined the constant 



^ ^^^^ 
where we recall that ATq = \m\{l — k)c^ is the freezing range. Note that the parameter H 
can be expressed in terms of the surface tension, H = '-//{IW). Then, we have 

A = IAToW/{2T), (60) 

where F is the Gibbs-Thomson constant of Eq. p6p . Therefore, up to numerical constants, A 
is the dimensionless ratio of interface thickness times freezing range and the Gibbs-Thomson 
constant. It is immediately clear that a variation of the interface thickness corresponds to 
a change in A. 



B. Non-variational formulations 



In spite of the theoretical appeal of a variational formulation, relaxing the requirement 
that both Eqs. ()31|) and ()32j) derive from a single functional F yields more flexibility. In 
particular, this extra freedom can then be used to cancel out spurious effects. 

a. Non-variational formulation without antitrapping current. In the last form proposed 
in the previous section, the interpolation function g{(j)) enters the model not only in the 
evolution equation for the phase field [Eq. (jSHI)], but also in that for the impurity, Eq. 
through the change of variable Eq. (j57|) . Whereas the condition ^'(±1) = is necessary 
in the equation of motion for (p to ensure that = ±1 are the equilibrium solutions for 
arbitrary u and T, no such condition is needed in the equation for the impurity. This 
suggests replacing g{(f)) in the definition of u Eq. (j37j) by another function h{(j)) which not 
necessarily satisfies /i'(±l) = 0, but still has the same limits h{±l) = ±1: 

" = '°( c?|i + t--(i-^M^)] )- 

Thus, the equilibrium properties derived in the last section are preserved; note, however, 
that the equilibrium concentration profile Co(0) is modified because h{(p) replaces ^(0) in 
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Eq. (f^ . yielding Co(0) = c°[l + — (1 — k)h{(j))]/2. In practice, this allows the simple choice 
h{(f)) = 0, for which the equilibrium concentration profile has the lowest possible gradients, 
and convergence of the simulations can hence be achieved for a coarser mesh isl ]. 

b. Non-variational formulation with antitrapping current. Albeit now h{(j)) and q{(j)) 
are completely free functions which purely need to interpolate from +1 to -1 and from to 
1 respectively, this does not yet provide enough freedom to cancel the three spurious effects 
mentioned in the introduction. To achieve this goal, we add an extra term in the model 
equations to specifically cancel one of them. The extra interpolation function contained in 
this new term provides the necessary third degree of freedom to make all three effects vanish. 

We specifically target the solute trapping effect. This occurs when solute atoms or 
molecules cannot escape the advancing solidification front fast enough to maintain local 
equilibrium at the interface. The characteristic interface velocities where solute trapping 
becomes important can be estimated by comparing the time of advance by one interface 
thickness, W/V, and the time it takes for the solute to diffuse through the interface, W'^/D. 
The result is ^ ~ D/W, and hence the critical speed depends on the interface thickness. 
Since we ultimately want to simulate solidification with diffuse interfaces that are orders of 
magnitude larger than the real solid-liquid interfaces, solute trapping sets in for much lower 
speeds than in reality. 

To eliminate this interface-thickness effect, we introduce a supplementary current in the 
equation for the solute concentration, the antitrapping current. Its purpose is to transport 
solute atoms from the solid to the liquid. Therefore, it has to fulfill a number of properties. 
First, it must be proportional to the speed of the interface, and hence to dt(f). Next, it must 
be directed from the solid to the liquid, that is, along the unit normal vector n, which in 
terms of the phase field can be expressed as (up to higher order corrections in the interface 
thickness) fi = — V0/|V0|. Furthermore, it must be proportional to the interface thickness 
W, and to the local concentration difference between solid and liquid. In contrast, we do not 
know a priori the profile of the current function through the interface. The time derivative 
of the phase field dt4> is sensibly different from zero only in the interface regions and induces 
a certain antitrapping current profile. Additional freedom may be gained by allowing for 
a shape function a(0) that must be appropriately chosen in order to obtain the correct 
thin-interface limit. 
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In summary, we write 



jat = aima - k)cy^n = -a{<P)W{l - A:)c?e"^^, (62) 

and the equation for the concentration becomes 

dc 



V(Dq{ct>)cVu-jat). (63) 



dt 

Note that the latter no longer derives from a functional F, even if such a functional is allowed 
to be different from that giving rise to the equation of motion for 0. 

c. Formulation in terms of dimensionless supers aturation and relation with pure sub- 
stance model. It turns out to be advantageous for the subsequent asymptotic analysis to 
make another change of variables in order to bring the equations in a form that is close to 
those analyzed in Refs. [isl . llS^ . To this end, we introduce the diffuse-interface extension 
f/(0) of the dimensionless supersaturation U in Eq. now defined in the whole system, 

- 1 

U^—. (64) 

Furthermore, we fix now the interpolation function g to be 
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A3 



9W = yU-^ + ^ . (65) 



define new interpolation functions 



, A + k - [l - k)h((j)) , , 

m = m — (66) 

^(0) = ^^(0)=(0-^ + f) (67) 

and transform the equation for c into one for U. Taking into account that T{z) = Tq + 
G{z — Vpt) and the temperature-dependent relaxation time t = tq[1 — {1 — k){z — V^t)//^] 
discussed later in section IV. C, the final set of equations is 



To 



T 



^ = W'V'<P + <P-<P'-X g'{<P) {U + '—^) , (68) 



- '-^m] f = V ■ D,i^) w + aim [1 + (1 - m 



ll + (l-.)f/]i^. ' (69) 



22 



where 



A = ^A. (70) 



With these choices, the phase field equation [Eq. (|68|) ] becomes identical to the one analyzed 
in Ref. One important advantage of this formulation is that the special case of a 

constant concentration jump can be recovered without any difficulty by setting k = 1, 
whereas in the formulation with the variable u, the limit k 1 has to be treated with some 
ca.e. Hence, the ,„odel of Ref. Q is e„nta,ned a. a speca. ease of E,s. 63 and 63. for 

k = l. 



IV. THIN-INTERFACE ANALYSIS 



A. Introductory remarks 

The goal of the matched asymptotic analysis is to relate the phase-field model [Eqs. 
()(i8|l and to a free-boundary problem. In particular, we would like to recover that 

of Eqs. ()26fOH|) . The principle is to choose the interface width much smaller than any 
physically relevant length scale. This difference in scale can be exploited for a perturbation 
expansion, in which the solution on the outer scale of the transport field is first assumed to 
be known. For a given point of the interface, this fixes the local velocity and curvature. The 
reaction of the diffuse interface to this "forcing" can then be calculated on the inner scale 
of the interface width, which yields a boundary condition for the diffusion field on the outer 
scale. The matching of both solutions then provides the link between "outer" (physical) and 
"inner" (phase-field) parameters. 

Two different perturbation schemes have been used. The "classic" one, developed by 
Langer, Caginalp, and others, uses the ratio of interface thickness and capillary length, e = 
W/do, as an expansion parameter. Later, Karma and Rappel remarked that the physically 
relevant length scales for the outer problem are not the capillary length, but rather the 
diffusion length D/V or a local radius of curvature p. Calculations performed with the 
expansion parameter p = WV/D for the symmetric model of solidification [Dg = Di, or 
q{(f)) = 1) yield, to first order in p, a new expression for the interface kinetic coefficient 
that contains a finite-interface thickness correction. This has allowed a tremendous gain in 
calculation power, since much larger W, including e ^ 1, can be used. It was also shown 
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that this correction can be obtained in a second-order expansion in e |13l . Il8| . 

Here, we will follow the classic scheme and present the asymptotic analysis for our model 
up to second order in e. While e is not necessarily small, this method yields all important 
correction terms at second order, while other schemes need to include some third order 
terms. The reasons for this, as well as the conditions of convergence of the expansion in e, 
can be better appreciated in the light of the formal results given below, and a discussion of 
these points is therefore deferred to Sec. IV. D. 

To perform the analysis, it is advisable to use a dimensionless version of the equations. 
We will use as unit length the capillary length do and as unit time dl/D. Without loss of 
generality, we set t = (which amounts to a shift of reference frame) such that the term Vpt 
drops out. Furthermore, we remark that from the definitions of Eqs. (I21|) . (j6(J|) . and (|7U|) . 
we obtain 

W 

do = ai— (71) 

with Oi = I/J, where J = (?(+!) — g{—l). For our choices of functions, / = 2-\/2/3 and 
J = 16/15, such that ai = 5^/2/8. Therefore, A can be eliminated of the equations in favor 
of aie. The result reads 

ae'dt(j) = e^V^ - f'{<P) - a^eg'{<f)) {U + vz) , (72) 
- d,U = V- [q{ct.) VU + ea(0) [1 + (1 - k)U] ^ ^ 

+ + (73) 

where we have introduced the dimensionless parameters u = do/lx and a = Dt/W"^, and 
defined the double-well function / = —0^/2 + 0^/4. We will assume that e is the only small 
parameter and consider all other parameters of 0(1). Note that u = do/l-r is a physical 
parameter that is typically small, but independent of the computational parameter e, and 
therefore, u = 0(1). The parameter a depends on the choice of r; we consider it to be of 
0(1) in order to avoid neglecting any important terms. Our conclusions remain valid if a is 
of order e or smaller. 

For comparison, we also adimensionalize the free-boundary problem we would like to 
recover, Eqs. using the above rescaling of space and time: 

dtU = V^f/ (liquid) , (74) 
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[l + {l-k)U] V, 



'n 



dnU\'^ (interface) , 



(75) 
(76) 



U 



K, — [3vn — yz (interface) , 



where k = d^K, and f„ = doVn/D are the dimensionless interface curvature and normal 
velocity, and (3 = (3D /do, the dimensionless kinetic coefficient. In the following, we will 
show how to recover this model as closely as possible by choosing specific forms for the 
functions g(0), h{(j)) and a(0). 

B. Matched asymptotic expansions 

We make a perturbation analysis in powers of e in the inner region. 



and similarly in the outer region, = 0o + e0i + . . U = Uo + eUi + . . .. In the outer region, 
Eqs. ()72|) and (fTSj) can be expanded in powers of e in a straightforward manner. Since we 
have (^'(±1) = 0, = ±1 are stable solutions for the phase-field equations to all orders in e 
for any value of U and z. Therefore, the outer solution for the phase field is simply a step 
function, and the field U obeys the diffusion equation to all orders. 



where we recall that g(l) = and q{—V) = 1 for the one-sided model. Also, note that the 
local equilibrium condition for the concentrations will be satisfied at all orders to which U 
is continuous across the interface. 

In the dimensionless equations, the Laplacian of the phase field comes with a prefactor 
e^, which leads to the two distinct constant solutions in the outer region on the two sides 
of the interfaces. In the inner region, the phase field varies smoothly. Equation (|72|) tells 
one that, for e ^ 0, this is only possible if such a variation takes place precisely on a scale 
of 0(e), which renders = 0(e~^) and invalidates the counting of orders used above. To 
compute the inner solution, we therefore must rescale the coordinate normal to the interface. 
We introduce the curvilinear coordinates in the reference frame of the interface r (signed 
distance to the level line = 0) and s (arclength along the interface), and define the rescaled 



U 



00 + e0i + + • • • , 

Uo + eUl + e'f/s + . . . , 



(77) 
(78) 



dtU = q{±l)V^U, 



(79) 
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coordinate rj = r/e. Standard formulas of differential geometry yield [see e.g. 

dt = -e'^Vndn + dt-vtds + 0{e), 

V2 = + e-'t^d, - ft:Vr, + d^, + 0(e), 

V ■ (gV) = e^'^d^iqd^) + e'^nqd^ - K^qr]dr, + d^igds) + 0(e), 

z = Zi + e{n- z)r], 

= n[l + 0(e2)l + 50(e), 

|V0| L ^ 

V ■ a = e~^drj{n ■ a) + ds{s ■ a) + kti ■ a + 0(e), 



where a is a vector function of the fields, Vn{vt) are the dimensionless normal (tangential) 
velocity of the interface, zi its dimensionless z position, and dt is the time derivative at fixed 
r and s. 

Since changes in the arclength s amount to a re-parametrization, we neglect terms in 
Vt without loss of generality. We will also neglect the operators dt. This amounts to the 
assumption that the interface follows adiabatically the changes in the forcing. For the 
phase field 0, this approximation is always justified, since this field has an approximately 
stationary kink shape moving with the interface (this will be explicitly checked by computing 
(j) at lowest order in e, which turns out to be a function of 7] only). For the diffusion field 
f/, it can be seen from Eq. ()76|) that dtU 7^ originates from variations with time of the 
interface curvature, velocity, and position. The variations of the latter occur generally on 
the slow time scale of solute redistribution transients, D /V^ , and are therefore negligibly 
small. The characteristic time scale for variations of the curvature and velocity is R/Vn, 
where i? = 1//C is the local radius of curvature, since this is the time the interface needs to 
move by once its radius of curvature. Therefore, the curvature and velocity terms in dtU are 
of order f„K(K + (3vn)- Since k and f„ themselves are small quantities, dtU is much smaller 
than other terms of order VnH which will appear in the calculation below, and can hence 
safely be dropped. 

We substitute the above expressions into Eqs. (|7^ and (f7H|) to obtain 

dl<P - f\<P) 
+ e [{avn + K)dr^(j) - ai5f'(0) {U + uzi)] 

'dlct> - ,^^vdr,<P - aiu{n ■ z)7]g'{<p)] = 0{e^), (80) 
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e-^dr,{qdr,U) 



+ t 



'l + k 1-k 



h{(f)) + nq 



+ Vndr, {a [1 + (1 - k)U] 9,0} - y [1 + (1 - k)U] S./lj 

+ e° [ds{qdsU) - K^mdr^U + aw„K [1 + (1 - k)U] 9,0} = 0(e) 



(81) 



and solve them order by order in e. The matching to the outer expansion is trivial for 
since the outer solution is just a step function. For f/, the matching conditions read 



lim 

r7=itoo 

lim 

»;=itoo 

lim 

»7=itoo 



f/o(r7, s) - f/o| 



0, 



U,{7],s)-m^{s)+T^drU^\^{s) 



0, 



0, 



(82) 



where |^ means that the outer field and its derivatives are evaluated at the interface, coming 
from either the + (liquid) or the — (solid) side. As a consequence, 

lim dnUo{r],s) = lim dlUi{r], s) = 0, 



lim drjUi{r],s), 

ri=±oo 



lim 

»7=itoo 



d,U2{v,s,t)-r]d^^Uo\^{s] 



(83) 



This matching will provide the boundary conditions on the interface for the outer concen 
tration. We now proceed to solve the inner equations order by order. 
Gibbs-Thomson relation. Equation (|8(J|) at order e°, 

d'jo - /'(0o) = 

yields, with the boundary conditions 0o —1 for r] —>■ +oo and 0o 
the matching to the outer solution, the zeroth order solution 

V 



1 for 7] 



(84) 

oo set by 



^o(^) = — tanh 



V2- 



In turn, Eq. (jHT|) at order e ^ becomes 

dr, {q{(t>o)driUa) = 0, 

which can be integrated once to yield q{(f)o)drfUo = Aq{s). Taking the t] 
according to Eq. (jH!I|) we find Ao{s) = 0, and therefore 

Uo = Uo{s). 



(85) 

(86) 
±oo limit 

(87) 
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To fix this constant, in turn, we consider Eq. (jHUj) at order e, 

£01 = aig'{(f)o) {Uo + vzi) - (aw„ + k)9^(/'o, (88) 

where £ = (9^ — f"{(f)o) is a hnear differential operator. Since the partial derivative with 
respect to r] of Eq. (jH^ is £9^0o = 0, (9^0o is an eigenfunction of C with eigenvalue zero. 
Therefore, the solvability condition for the existence of a nontrivial solution (pi reads 

ai (Uo + uzi) J + (avn + k)I = 0, (89) 

where J = g'i(po)dri<Pod'n = gi+l) - gi-l) and / = {dr^cp^f drj. Since / and J are 
the same constants that have been used to define ai = J/J, we obtain 

Uq = —UZi — aVn — K, (90) 

which is identical to the Gibbs- Thomson condition of the free boundary problem, Eq. ()76p . 
with P = Po = a. 

This is the "classic" result for the kinetic coefficient in the sharp-interface limit. To 
obtain the thin-interface correction, we repeat the same procedure at next order. Thanks 
to Eq. (j87|) we can drop the terms in 9^?7o arising in Eq. (j81|) at order e~^, to obtain 

dr, [q{<Po)d,Ui] = -Vr^dr, {a(0o) [1 + (1 - k)Uo] d.M + y [1 + (1 - k)Uo] dr,h{<Po), (91) 

and integrate it once with respect to rj to yield 

g(</.o)a,C/i =Vn[l + il- k)Uo] [/i(0o)/2 - a(0o)5^0o] + ^i(s), (92) 

where Ai{s) is an integration constant. The latter can be fixed by considering the limit 
7] ^ —oo. In fact, the left-hand side represents the diffusion current, which vanishes inside 
the bulk solid. Since the antitrapping current must also vanish and h{l) = 1, we find 
Ai{s) = — (f„/2) [1 + (1 — k)UQ]. Substituting it back into Eq. (jHH and integrating the 
latter once more between and t], we find 

Ui = Ui + ^[l + il- k)Uo] PpIMO] d^, (93) 

Z JO 

where Ui is the value of Ui at the interface {rj = 0), and 

/ , N h{(Po) - 1 - 2a(0o)(9^0o 

= m ■ 
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The profile Ui therefore depends on the choice of the functions q{(f)), h{(j)) and a(0). Note 
that both the denominator and the numerator tend to zero when 7] —>■ — oo. It is important 
here to remark that we need to require p ^ in this hmit, since otherwise Ui diverges, which 
makes a matching to the outer solution {U is constant in the solid) impossible. In fact, this 
property makes the standard asymptotic expansion inconsistent. A careful analysis, carried 
out in the appendix, shows that in this case a term of order plogp (with p = WV/D the 
interface Peclet number) appears in the interface kinetics, which makes the convergence of 
the model to the sharp-interface limit very slow. This term appears, for example, in the 
standard formulation of the one-sided model that has been widely used [3, 12|- In order to 



avoid this phenomenon, we will require in the following p(0) for — >■ 1, that is, the 
numerator must vanish more rapidly than the denominator. 

Under this condition, we may fix the constant f/i by considering Eq. ()80|) at order e^. 



C-(t)2 = — ^^0i-(at;„ + K)9^0i + ai5('(0o)f^i + fl'"(0o)0i«i(t^o + '^^i) 

+ n'^'ndncpo + g'{(f)o)aiu{n ■ z)r] (95) 

where we have used ds4>o = 0. In this expression appears the first order correction to the 
phase-field, 0i, which is the solution of the differential equation obtained by substitution of 
Eq. dHl into Eq. (jHHI), 

£01 = -{avn + k) (ai5('(0o) + 5r,0o) , (96) 

with the boundary conditions (f)i{ri — > ±oo) = imposed by the matching to the outer 
solution. Clearly, 0i equals aVn + k times a function only of rj, so that, when substituted 
into Eq. it would yield {aVn + contributions to tJi. There are essentially two ways 
to avoid this problem. The first would be to choose g such that 

g'{(f)o) = -a,0o/ai, (97) 

which makes 0i vanish. For our standard quartic double-well potential which yields 
dri(po = (1 ~ 0o)/v^, the corresponding g function is a third-order polynomial that has 
been widely used. However, we have chosen here a different function, and many calculations 
have also been performed with yet other interpolation functions, such that this condition 
is too restrictive. The second way out is to use the symmetry properties of the involved 
functions. For any symmetric double- well function (that is, /(— 0) = /(0)), the equilibrium 

29 



profile is odd in 77, (pQ^—r]) = —cpQ^r]), and its derivative is even. If g is cliosen to be odd 
in 0, (?(— 0) = —g{(j)), tlien g'{4'o) is also even in rj. Therefore, the entire right-hand side of 
Eq. (jnSl) is even. Since C is also an even operator, 0i must be even, and its derivative 
odd. Given that the solvability condition is obtained by multiplying the right-hand side of 
Eq. by 9^0O; an even function, and integrating from —00 to +00, the contribution of 
all odd terms vanishes. The only remaining is the one that contains Ui, and the solvability 
condition reads 

^[l + {l-k)Uo]K~JU^ = (98) 
where we have expressed Ui according to Eq. (jHS)), and 

dvd,<Po9'{<p') / p{MO)dt (99) 

-00 Jo 

To obtain the desired result, namely an expression for Ui, let us first remark that in 
the limit 77 00, Eq. (jHS} yields 9^f/i = —Vn [1 + (1 — k)Uo], which is just the Stefan 
condition at lowest order. Using the matching conditions lim^^-too dr,Ui = dr Uq , and 
= lim^^-too Ui{ri) — ridrUol"^, we obtain 

U,\^ = -vJt. (100) 
TF^ -\- K 

f^f = -[! + (!- km J , (101) 
[p{<Po)-p{<P^)]dv. (102) 

Jo 

Note that U will be continuous across the interface up to 0(e) if and only if = = F 
(and hence Pi = /3f ). Since U = Uq + eUi, the total kinetic coefficient is 

~ ~ ~ K 4- JF^ 

f5^ = f3o + ef3t = a - e [1 + (1 - k)Uo] . (103) 

The implications of this finding will be discussed below. 

Mass conservation. As already mentioned before, Eq. (j^JHj) together with the match- 
ing condition (j83|l yields drUo\~ = and drUo\'^ = — f„ 1 + (1 — k)Uo , which is just the 
Stefan condition at lowest order. In order to evaluate eventual corrections, we proceed 
by calculating the normal gradients at order e using the matching condition for drUi\^ in 
Eqs. (I83|) . The quantity O^UqI^ can be evaluated by remarking that the outer problem sat- 
isfies a simple diffusion equation in a moving curvilinear coordinate system, and therefore 
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[drr + {vn + i^)dr + dss]UQ = 0, such that d^Uol^ = -[{vn + K)dr + dss]Uo\^- To obtain dnU2{r]), 
Eq. (|HT|) is evaluated at 0(e°) and integrated once from to r], 

+ kJ^' d^q{(j)o)d^Ui+VnJ^' - ^-^/i(0o)) d^Ui 

+ {a'(0o)0i [1 + (1 - A;)f/o] + a(0o)(l - k)Ui} 9,0o 

+ t;„a(0o) [1 + (1 - A:)f/o] 9^01 + VnK [1 + (1 - A;)f/o] T rf^ a(0o)95</)o 

Jo 

- ^(1 - A;) r rfe [719^/1(00) - ^ [1 + (1 - A;)?7o] /i'(0o)0i 

+ 9,,t/o r d^q{<Po) = A2{s), (104) 

JO 

where ^2(5) is an integration constant and we have taken into account that 9^f/o = ds(f)o = 0. 
Fortunately, we can drop many terms of this long equation because we are only interested in 
the limits 77 ±00. In this limit, 0i and (9^0o are exponentially small, such that all terms 
containing them can be dropped, except when they appear under an integral. The third 
term on the left-hand side of Eq. p()4p can be rewritten using Eqs. and ()94|1 as 

(3) = ^ [1 + (1 _ k)Uo] di [/i(0o) - 1 - 2a(0o)5c0o] , (105) 

and it can be seen that the part proportional to a(</'o) cancels out with the seventh term 
on the left-hand side. The remaining piece can be rewritten, using the Stefan condition to 
lowest order, as 

hm [(3) + (7)] = Kr,drU,\^ + ^ [1 + (1 _ k)U,] T di [/i(0o) - KtI)] ■ (106) 

Next, the remaining terms that contain h can be grouped and integrated to yield 

(4) + (8) = - ^^M0o)) UM- (107) 

Using the matching condition for f/i, lim^^-too t^i(^) = +V'^rUo\^, and the fact that 
lim^^_oo q{(t>o)driU2 = 0, we can obtain the constant A2 from the limit r] —00 of the entire 

Eq. (uni, 

= ^ [1 + (1 - k)Uo] / d7] [/i(0o) - 1] + Vr,kU{ + dM / dr]q{<f)o) . (108) 
2 JO JO 

Next, lim^^oo Q'(0o)t^r?f^2 is evaluated using the above result for A2. Finally, with the help of 
the matching condition and the expression for 9^t/o|^, we obtain 

drUi\^ = [1 + (1 - k)U,] {H+ - H-) - dMQ^ - Q-) 

- Vn{l - k)U+ - Vnk{U+ - U^) (109) 
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with 



±00 



±00 



(110) 

dv[qiMv))-qi<t^^)]- (m) 

The first two terms on the right-hand side of Eq. ()1()9|) are the announced finite interface 



thickness effects associated with interface stretching and surface diffusion; the third is the 
expected first-order term that appears on the left-hand side of the Stefan condition, Eq. (f73j): 
finally, the last one is a correction associated with a jump of U through the interface. In 
total, the mass conservation condition for the outer fields up to first order reads (recall that 

Uo = Uo) 



1 + {1 - k){Uo + eU, 



-driUo + eU,) + e 

+ dssUo{Q^-Q 
vlk 



+ 



l + {l-k)Uo (F+-F-) . 



;ii2) 



C. Discussion 



Physical interpretation of the corrections: There are three corrections in e to the classic 
free-boundary problem. The term proportional to — Q~ describes the response of the 
interface to lateral concentration gradients, caused by variations of the curvature or the 
growth speed along the interface. For a diffuse interface, the resulting mass flow is smaller 
than in the bulk liquid on the liquid side, but larger than in the bulk solid on the solid side. 
If the two effects do not exactly compensate, a surface diffusion term needs to be included 
in the Stefan condition. The condition to make this correction vanish is = Q~, which 
can be shown to be exactly the same as Eq. in the introduction by taking into account 
that g(0) = g(0)co(0)/Q. 

Next, the term proportional to — H arises from the source term in the U equation. 
If a positively curved interface moves forward, the liquid side of the interface is slightly 
longer than the solid side. Therefore, the source term on the liquid side is active over a 
larger area than the one on the solid side, and the integral of the source strength multiplied 
by the area over which it is active is precisely given by the difference — H~ . If this 
quantity is non- vanishing, the interface acquires a "net impurity content" , that is, a source 
term appears in the mass conservation condition when the length of the interface changes, 
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which is precisely the case if the product f„K is non-zero. This is the interface stretching 
correction, which vanishes when = H~ . In terms of the concentration, this condition is 
identical to Eq. 

Finally, the last correction involves a macroscopic discontinuity in U that is proportional 
to the velocity f„ and to — F~ , and that appears in the boundary conditions at the 
interface and in the Stefan condition, Eq. ()112|) . This is the solute trapping term: since 
the concentrations on both sides of the interface vary with velocity, they do not satisfy the 
partition relation Cg = kci out of equilibrium, or, in other words, the solute rejection is 
velocity-dependent. Since U can be assimilated to a chemical potential, its jump can be 
interpreted as resulting from a finite interface mobility that leads to interface dissipation. 
Note that both analogies are limited: whereas a "physical" dissipation is necessarily positive, 
the difference F^ — F~ here can have either sign, depending on the choice of the interpolation 
functions. Without the antitrapping current (a(0) = 0), the condition F^ = F~ that makes 
this correction vanish is identical to Eq. ^ in the introduction. 

Choice of functions: In order to make all three corrections cited above vanish, we need 
to satisfy simultaneously three conditions, namely, 

F+ = F- H+ = H g+ = Q-. (113) 

For fixed double- well and tilting functions / and g, we have at our disposal three interpo- 
lation functions: the diffusivity the source function h{(j)) and the antitrapping current 
profile a(0). The new element here is the antitrapping current. If it is absent, only two 
interpolation functions are available. It is then, of course, easy to satisfy two out of the 
three conditions. For example, choosing h odd in and q{(j)) = 1 — q{—(f)), respectively, 
will automatically satisfy the interface stretching and surface diffusion conditions. However, 
as already discussed in the introduction and also by Almgren for a thermal model jl^, all 
three of them can be satisfied only for a weak contrast in the bulk diffusivities, which of 
course excludes the one-sided case of interest here. The problem is that, in order to satisfy 
the integral conditions shown above, the interpolation functions need to be non-monotonous 
or even to change sign, which leads to strong higher-order correction terms or even to the 
emergence of singularities. 

It is interesting to note here why the corrections to the Stefan condition, namely interface 
stretching and surface diffusion, which were not computed in Ref. [l3|, vanish for the 
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symmetric model of solidification, g(0) = 1, a(0) = 0. Obviously, surface diffusion does just 
not arise for a constant q{(f)). But, furthermore, reduces to h{(j)) — 1, and therefore the 
two conditions = F~ and = H~ become identical, such that the "miraculous" choice 
of M.) odd i„ , which e„su..ed ^ F i„ Ref. Q also canoe, .he h.erface ..etching 
correction. 

The more involved one-sided case is cured with the help of the antitrapping current, 
which offers an additional degree of freedom to satisfy the third condition. The only place 
where the function a(0) appears in the final results of the matched asymptotics is in the 
first-order concentration profile Ui, and more precisely in the function = [h{(j)) — 1 — 
2a(0)(9^0o]/Q'(0)- A suitable choice for the function a(0o) can be obtained by a simple 
analogy with the symmetric model of solidification. For the standard choices a(0) = and 
h{(f)) = (f), we have p{(j)) =0 — 1. The same function p{(j)) can be recovered in the one-sided 
case if we choose 

g(0) = (l-0)/2, (114) 
h{<P) = 0, (115) 
a(0) = (116) 

since we can exploit the fact that 9,,0o = — (l/v^)(l — 0o)- Then, all the solvability integrals 
are identical to those calculated in Ref. (isf for the symmetric model in the isothermal 
variational formulation. 

Essentially, this "trick" solves the problem because it makes the two conditions = 
F~ and = H~ identical, as for the symmetrical model. The same strategy can be 
applied to obtain other possible phase-field formulations. For any "source function" /;,(</)) 
and diffusivity g(0), the equivalence to the analogous symmetric model can be obtained by 
requiring p(0) = /i(0) — 1, which yields 

V2(02-l) ^ ' 

For example, the function f/i of the symmetric model in the variational formulation of 
Ref. Q|, which uses h{(j)) = g{(f)) = 15(0 — 20"^/3 + 0^/5)/8, can be recovered for g(0) = 
(l-^)/2by 

aW = K5^1^M%±lM±i. (US) 

2v2 

Since this model is known to be less efficient, we have not investigated further this alternative. 
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Kinetic coefficient: For low-speed solidification, kinetic effects are usually negligibly 
small, and therefore we want to make the kinetic coefficient vanish. This is possible because 
it consists of two contributions of opposite signs. Converting Eq. (|1(J3|) back to dimensional 
units, we find (in the following, we will assume = = F) 



T I XVV^ 1 

P = ^^xwY ~ + ~ ' ^^^^^ 

a. ^ (120) 

For k = 1 (constant concentration jump), this is identical to the expression for the symmetric 
model, and /? = can be achieved by choosing A = {jD) / {a2W'^). For A; 7^ 1, the kinetic 
coefficient depends on t/o, the average value of U in the diffuse interface. The physical 
meaning of this dependence can be understood as follows. The second term in the expression 
for (3 arises from the additional driving force supplied to the interface by the redistribution 
of solute inside the diffuse interface. For 7^ 1, the amplitude of this redistribution depends 
on the local state of the interface, since the concentration jump depends on temperature, 
curvature, and kinetics. To see this, recall that Uq = —Zi/lx — d^K — /3oKi, where /3o = 
aiT/{\W), according to the dimensional version of Eq. ()9()|1 . and furthermore that q/c^ = 
1 + (1 — fc)f/, and the concentration jump from solid to liquid is Ci{l — k). 

As a consequence, the interface kinetics depends on the local geometry and velocity of the 
interface, and it is not possible to make /3 completely vanish by the same choice as before. 
Among the correction terms, c/q/C and (3oVn are usually small, but no general statement can 
be made about the magnitude of zi/It- Two strategies are possible to tackle this problem. 
The first is to choose a temperature-dependent phase-field relaxation time, 

T = TQ[l-{l-k)z/lT]. (121) 

This does not change the asymptotic analysis for the equation since the z-dependent part 
does not contribute to the solvability conditions. It is sufficient to replace r by Eq. ()121|) in 
Eq. ()119|) . With the usual choice A = {tqD) / {a2W'^) , the residual kinetic coefficient is 

(3 = (3o{l-k){doK, + (3oVn), (122) 

with I3q = aiTo/{XW). The temperature-dependence is eliminated, but curvature and veloc- 
ity corrections to (3 remain. 
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The second strategy is to introduce a ^/-dependent phase-field relaxation time, 

r = To [1 + (1 - k)U] . (123) 

The idea is to make both terms of Eq. ()103j) contain the same prefactor [1 + (1 — k)Uo\ such 
that the compensation of the two terms is independent of Uq. This time, the solvability 
conditions for 0i and 02 are modified. The former yields a the new expression for Uq 

1 + avn[l — k) 



Equation ()98jl that yields Ui becomes 

ai|^[l + (l-A;)f/o]ir- Jf/i} 

- avn{l - A;) [1 + (1 - k)Uo] |/f/i - ^ [1 + (1 - k)U,] K'^ = 0, (125) 

where the new solvabihty integral, 

/OO pT} 
dvidM' p{MO)d^, (126) 
-OO Jo 

equals K' = 0.1869 for the choice of interpolation functions given above. A straightforward 
calculation yields 

,7 _^"M^n 1-^( 1 + - fc) [1 + (1 - k)Uo] [K'J/{KI)] ] 

An expansion of this result in Vn shows that the leading order prefactors of the two terms in f3 
originating from Uo and Ui are indeed the same. Furthermore, it can be seen that all higher 
order corrections are proportional to aVn{l — k) = /5oV^(l — k) = [aiT^ / {W X)]Vn{l — k). As 
long as this quantity is much smaller than unity, the resulting residual kinetics should be 
small. 



D. Limits of validity and expansion parameters 

In the numerical calculations presented below, we obtain converged quantitative results 
for values of e = W/do much larger than unity, even though we have used e as a small 
expansion parameter in the thin interface analysis. This raises the question: what is the 
domain of validity of this expansion ? A rigorous answer to this question would in principle 
require to carry out the expansion in e at one more order to determine when the additional 
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corrections to the boundary conditions are negligible for a given set of growth conditions. 
This represents a formidable analytical task that is beyond the scope of this work. We can, 
however, use dimensional arguments to place bounds on the validity of the thin interface 
analysis. We shall conclude from the foregoing arguments that e need not be small for this 
analysis to be valid, consistent with the numerical findings; W only needs to be smaller than 
a characteristic length do, which depends generally on the growth conditions. 

The expansion defined by Eqs. ()77p and ()78|) assumes that e is small and that the func- 
tions (j)n and Un are of order unity. The magnitudes of the functions [/„, however, are not 
known without specifying the outer problem. For typical growth conditions, the variation of 
concentration along the interface due to capillarity and interface kinetics is small. In partic- 
ular, the velocity-dependent form of the Gibbs-Thomson condition implies that \U + vz\ <^ 1 
in the diffuse interface region as long as k + (3vn <^ 1, and that therefore the right hand side 
of Eq. (f7^ contains small terms other than e. To define a diffusion field that is of order 
unity in the interface region, consistent with the choice of e as a small expansion parameter, 
one would need to rescale the combination U + vz inside the interface by some character- 
istic mean interfacial value of the diffusion field, f/, which depends on the outer solution. 
This procedure, however, does not change the results of the asymptotic analysis because it 
amounts to a simple change of variable. For convenience, we have therefore opted to keep 
the expansion parameter e, which is independent of the outer solution. It is clear from the 
above arguments, however, that this expansion is valid as long as 

ef7 < 1. (128) 

Since IJ is typically small, e need not be small for the expansion to be valid. 

To make Eq. ()128j) more transparent, it is useful to re-express this constraint in terms of 
the interface velocity Vn and the local radius of curvature R. Up to coefficients of order unity, 
which we do not consider, and assuming that the velocity is positive, it follows dimensionally 
that t7 ~ |f/ + z/2;| ~ do/ R + [3Vn + WD/Vn, where do/R and pVn are capillary and kinetic 
corrections originating from the velocity-dependent form of the Gibbs-Thomson condition, 
and WD/Vn originates from solute diffusion in the diffuse interface region. The product 
eU is therefore of order W/R + (3VnW/do + W'^Vn/ (doD). In terms of e, the dimensionless 
kinetic coefficient (3, and the Peclet number p, Eq. ()128j) can be rewritten as 

W/R + p{(3 + e) <^1. (129) 
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The same estimation can be obtained directly from the expressions for Uo and Ui calculated 
above. Convergence is hence limited by two independent conditions, linked to the local 
curvature and velocity, respectively. The first condition, W/R <^ 1, states that the interface 
thickness must be much smaller than the local radius of curvature. The interpretation of the 
second condition, p{(3 + e) ^ 1, depends on the physical value of the kinetic coefficient to 
be simulated. In the present work, we focus on the limit of vanishing kinetic effects relevant 
for small growth velocity (/3 = 0), which is achieved by setting r ~ XW'^/D. Therefore, the 
limiting condition is pe <^ 1, which can also be rewritten as rVn/W -C 1. In practice, we 
found that the convergence starts to break down for rVn/W ~ 0.2 or W/R ~ 0.2, although 
occasionally slightly larger values of rVn/W could be used. 

Defining the diffusion length I = D/Vn, Eq. (jl29j) can also be rewritten in the form 
W/ic 1, where ic = do/ {do/ R + (3Vn + W/l). This shows that the true small parameter 
elJ can always be expressed as the ratio of W and some characteristic scale scale ^c which 
is much larger than d^ and which depends on experimental growth conditions. 

Finally, it is in principle possible to use the interface Peclet number p = WV/ D as a small 
expansion parameter in the thin-interface analysis, as for the solidification of pure melts with 
symmetrical diffusion However, this choice is not optimal for the case of asymmetrical 
diffusion considered here for technical reasons. In particular, the interface stretching and 
surface diffusion terms appear at second order and third order, respectively, in an expansion 
in p. In contrast, they appear both at second order in the e expansion. Therefore, the latter 
is preferable for clarity of exposition, with the caveat that it is necessary to consider the 
outer region to obtain the true condition of validity of this expansion expressed by Eqs. 
(inHll-dnni), or equlvalently by the condition W/£c < 1. 



E. Anisotropy 



To include anisotropy, it is sufficient to proceed in the standard manner, that is, make 
W and T orientation-dependent, as in Refs. 13, 

\4' 



W{n) = Wa,{n) = W{1 - Se^) 



l-3e4 



iv0r 



T(n) = roa^(n) 



(130) 
(131) 
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Here, it is understood that tq might be replaced by its temperature- or [/-dependent version. 
As a consequence, the standard resuh for the anisotropic capillary length is recovered. 
For the interface kinetics, the orientation dependence appears together with tq in all the 
above results. Finally, note that the interface thickness also appears as a prefactor in 
the antitrapping current. However, since the anisotropy of W itself is small (recall that 
the anisotropy of the capillary length is 15 times larger that the one of W for fourfold 
symmetry), only a small error will be made if the actual orientation-dependent interface 
width is replaced by its mean value in this term. 

V. NUMERICAL TESTS 

We have simulated the phase-field model of the directional solidification of a dilute binary 
alloy defined by the anisotropic version of Eqs. (|68l) and (|69|) for parameters corresponding 
to the impure succinonitrile (SCN) alloy of Ref. (22j. The alloy parameters together with 
the values of the pulling speed and the temperature gradient are listed in Table HI The 
chosen pulling speed is ten times the value for the onset of the MuUins-Sekerka instability. 
For these parameters, the capillary length is several orders of magnitude smaller than the 
thermal length or the diffusion length. Since typical cell widths are ~ 100 fim or ~ lO^do 
and computations are only feasible if one cell width ~ 10^ grid points, we are forced to use 
values of W much larger than do, typically W/d^ — 10 to 100. We will see that, with the 
present phase field model, it is possible to obtain well converged results even with such large 
W/do ratios. 

To choose the phase field model parameters, we first note that the ratio of the capillary 
and thermal lengths, v = d^/lx = 4 x 10^^, and the dimensionless pulling speed Vp = 
Vpdo/D = 4.16 X 10~^ completely specify the interface evolution in the sharp- interface 
equations. This can be seen by scaling length and time in these equations by do and d^/D, 
respectively. In the phase-field model, we have the additional length W and converged results 
should be independent of the ratio e = W/do- Note that for anisotropic surface tension 
Vr(n) = Wasin) with as{n) given by Eq. ()130p . In a given simulation, we fix e = W/do 
and hence A = aie from Eq. (ffT| . Furthermore, we use a temperature- and orientation- 
dependent relaxation time r as specified in the previous section together with the relation 
To = a2\W^ / D , which makes the interface kinetic coefficient vanish for all temperatures and 
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m Coo (shift in melting temperature) 


2 K 


D (diffusion coefficient) 


10-9 mVs 


r (Gibbs-Thompson coefficient) 


6.48 X 10-*^ Km 


Vp (pulling speed) 


32 fim/s 


G (thermal gradient) 


140 K/cm 


do (capillary length) 


1.3 X 10-2 fim 


It (thermal length) 


3.33 X 10^ ^m 


Id (diffusion length) 


60 iim 


k (partition coefficient) 


0.3 



TABLE I: Parameters for the impure succinonitrile (SON) alloy system of Ref. |27| used in the 
phase-field simulations and corresponding characteristic length scales for directional solidification. 
The anisotropy of the interfacial free energy is taken to be £4 = 0.007 (0.7% anisotropy). 

orientations, and we scale lengths by W and time by tq in the phase-field equations. The 
scaled phase-field equations then only depend on e through the dimensionless parameters 
D = Dtq/W"^ = aia2e, Vp = VpTo/W = Vpaia2e'^, and It = It/W = l/(ez/). Writing out 
explicitly all the interpolation functions, and taking into account the contributions of the 
anisotropic W{n) in the functional derivative, the equations read 
z — V„t 



{i-kY- 



ii 



a An) — = V 



dt 



a,(n)'V0 



+ 5x 

+ 0- 



V0 



da An) 



'l + k l-k 



'dt 





[1 + (1 - k)U] 



(133) 



where x and z are in units of W and t is in unit of tq. Simulations are repeated with 
different values of e to study the convergence. The equations are discretized on a square 
lattice; some details are given in the appendix. We have used a grid spacing Ax/VT = 0.8 
in most of the simulations, but also used a finer resolution /S.x/W = 0.4 to study the effect 
of the discretization. For the time evolution, we have used an explicit Euler scheme with a 
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time step chosen below the threshold of numerical instability for the diffusion equation in 
two dimensions, At < (Ax)^/(4D). 



A. Stability spectrum 




FIG. 2: Comparison between the linear stability spectrum of a planar steady-state interface com- 



puted with the phase-field mode^ 
and the Mullins-Sekerka theory 



for different interface thicknesses (dot-dashed and dotted lines) 



2J] (solid line). Here, x is the growth rate of a sinusoidal per- 



turbation of wave number Q, and / = 2D/Vp is the diffusion length. The parameters are for an 
impure SCN alloy system described in the text with Vp = 32 //m/sec and G = 140 K/cm. 



We have numerically calculated the stability spectrum of a planar steady-state interface. 
To this end, the system was initialized with a planar interface at its steady-state position. 
The concentration in the liquid was set to the exponential steady-state solution of the 
free boundary problem. A small sinusoidal perturbation of amplitude A <^ W and wave 
number Q was then applied, and its time evolution was followed by extracting successive 
interface positions. It follows an exponential increase or decay, and the growth rate x{Q) 
was extracted by a fit of the perturbation amplitude versus time. 

In Fig. El we compare the results from the numerical simulations to the analytical solution 
for the Mullins-Sekerka stability spectrum of the free-boundary problem of Eqs. (f7^ - ()76p . 
The convergence is better for smaller wave numbers, which is perfectly reasonable since the 
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FIG. 3: Convergence of the growth rate x{Q) as a function of W/do for: a) Ql = 10.5, and b) 
Ql = 87.3. The dotted lines are the predictions of the Mullins-Sekerka analysis. 

ratio of perturbation wavelength to interface thickness scales with Q. For W/d^ = 9.025, 
the phase field model gives a good agreement for almost the whole range of wave numbers, 
including the maximum, which is the most important part of the spectrum. In Fig. IHl we 
plot the growth rate x{Q) of two selected modes versus the ratio W/do, which shows a fast 
convergence. For Ax/W = 0.4, the results are fully converged to the theoretical value for 
W/df) = 4.51 even for the mode with high wavenumber. It can also be seen that the larger 
grid spacing of Ax/W = 0.8 introduces slight corrections that are due to the lattice pinning 
effect (see Ref. |l8|). 
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FIG. 4: Convergence of the shape of steady-state deep cells as a function of interface thickness. 
Lengths are scaled by the cell spacing A = 22.5 fj,m, Vp = 32 ^m/sec and G = 140 K/cm. 
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B. Cell shapes 

To asses the convergence of the models in the nonhnear regime, we have computed shapes 
of steady-state cells for various values of W/d^. The simulation box contains half of a cell, 
with no-flux boundary conditions along the cell center and the groove. We have considered 
narrow cells of spacing A = 1732.8 x do = 22.5 fim, since we want to test the convergence 
of the model for small tip radii; in an extended system, these cells would be unstable to a 
cell-elimination instability that leads to a doubling of the cell spacing. As initial condition 
we set Q = Coo, Cs = kci (which, with the definition of Eq. and using c° = Coo/k, 
corresponds to U = — 1 in the whole system), and add a small sinusoidal perturbation to 
the interface, with a wavelength equal to the cell spacing and its maximum located on the 
boundary. After a transient where the interface recoils, it reaches steady state in the form 
of a half cell. 

The resulting shapes are shown in Fig. EJ For the cell shapes the convergence is faster 
than for the growth rate, and already simulations with W/do — 50 are well converged. To 
show more clearly the difference in the speed of convergence, we plot in Fig. El the tip radius 
p/df), the tip undercooling f2 = 1 — Zup/lr (where z = corresponds to the position of the 
steady-state interface), and the solute concentration in the solid in the center of the cell. 
For the latter, we compare the values that are directly obtained from the simulations (that 
is, the value of the field U in the center of the cell) to the value expected from the Gibbs- 
Thomson condition and partition relation at the interface, cl^ = k[k + {1 — k){Q — d^/p)], 
where the values of Q and p are obtained from the numerical results (Figs. Et-b.) Again, 
all the quantities are well converged for W/do ~ 50, and even for the ratio W/d^ = 72.2 
(corresponding to p/W ~ 4), the error in the tip radius for the phase field model is only 
about 15%, while the equilibrium solute concentration condition at the interface is satisfied 
within an error of about 1%. The error in this latter condition is small, even for the largest 
values of W/do used (Fig. Efc). Since microsegregation is important for metallurgy, the 
precise calculation of the solute concentration in the solid is an important new feature of 
the present model. 
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FIG. 5: Convergence as a function of interface thickness of various quantities associated with 
steady-state cell shapes: a) tip radius p, b) dimensionless tip undercooling and c) solid con- 
centration in the center of the cell. The diamonds (dashed line) in c) correspond to the values 
calculated from the Gibbs-Thompson condition using the tip radius of the phase-field shape. The 
inset shows the relative error of the phase-field results with respect to the Gibbs-Thomson predic- 
tion. 

VI. CONCLUSIONS AND PERSPECTIVES 

We have presented a detailed asymptotic analysis of the phase-field model for alloy solid- 
ification that was introduced in Ref. 19|, and we have simulated directional solidification of 
a dilute binary alloy. We have found a very good quantitative agreement with the Mullins- 
Sekerka stability spectrum of a planar interface for typical experimental control parame- 
ters. For solidification cells, we found that the solute concentration inside the solid agrees 
self-consistently with the prediction of the Gibbs-Thomson condition, in contrast to earlier 
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models where the microsegregation was only qualitatively reproduced jl7|. 

This advance relies on a solution of the complete problem of canceling all relevant thin- 
interface corrections to the original free boundary problem. This opens the way for quanti- 
tative comparisons between experiments and simulations both in two and three dimensions, 
with the concomitant possibility of testing the theories and concepts used to interpret mi- 
crostructural pattern formation, as was previously done for dendritic solidification. 

The present work can be extended along several lines. For example, it has been demon- 
strated that the concept of the antitrapping current can be generalized to two-phase so- 
lidification, which makes it possible to study eutectic or peritectic composite growth with 
excellent precision [2^. Also, the present one-sided model can be combined with a symmet- 
ric thermal model to yield a quantitative thermosolutal model of solidification 2^ . A small 
solute diffusivity in the solid can also be introduced without appreciable modifications of 
the present analysis. Finally, the antitrapping current, which was used here to restore the 
equilibrium partition relation, can also be used to obtain a non-vanishing, specified trap- 
ping. This is especially important to extend this model to the whole range of solidification 
velocity relevant for experiments. In addition, the present model should be applicable to 
model Hele-Shaw flows when the viscosity of one fluid is much smaller than that of the other. 

From a broader perspective, this progress revives the hope of using the phase-fleld method 
as an efficient and fully predictive tool for other free boundary and interface growth prob- 
lems where the dynamics of the two media are not necessarily symmetric, even outside the 
framework of systems described by a Lyapounov functional. A key element of this progress 
is the use of non-variational terms which provide additional freedom to obtain the correct 
mapping between a diffuse interface model and a desired free-boundary problem, such as 
the antitrapping current here, and other terms in other contexts 2^. It is important to 
emphasize that the interface is spatially diffuse and all interpolation functions are smooth in 
the present phase-fleld model. Hence, this model remains simple to implement numeri callv 
in comparison to other methods that combine sharp and diffuse interface ingredients Isfll . 

Let us conclude with a few remarks on the formulation of the model itself. The ther- 
modynamic derivation presented here, which is an alternative to previous expositions of the 
same model Q], establishes new connections to other phase-fleld models of alloy solidiflca- 
tion. As mentioned before, early phase-fleld models of alloy solidiflcation were plagued by 
a dependence of the surface tension on the interface thickness that arose from the coupling 
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between the phase-field and concentration equations 0, 0] • This problem was solved later 
by the introduction of two separate concentration fields, one for the solid and one for the 
liquid, and by interpreting the interface as a mixture of two phases 2- The requirement of 
local equilibrium between the two phases then allows to eliminate one of the concentration 



fields 



The resulting model has a surface tension that is independent of the interface 



thickness and can be used for arbitrary phase diagrams; however, some thin-interface effects 
remain, in particular surface diffusion 

In our derivation, we have succeeded in constructing a quantitative model for an ideal 
dilute binary alloy with a single concentration field, but two different interpolation functions 
of the phase field for entropy and internal energy density. This is appealing from a ther- 
modynamic viewpoint, since it maintains the interpretation of the concentration as a local 
quantity rather than a two-phase mixture. An interesting task would be to generalize this 
approach to arbitrary phase diagrams and multi-phase solidification. 
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APPENDIX A: ANOMALOUS INTERFACE KINETICS 

We give here a more detailed discussion of the interface kinetics in a phase-field model 
without antitrapping current, and with source and diffusion functions given by = 
and q{(f)) = (1 — 0)/2. We will see that in this model, the interface kinetics has logarithmic 
corrections. This occurs whenever in the limit (p 1 the ratio {h — l)/q does not vanish 
(i.e., remains finite or diverges). Note that, in physical terms, the two functions describe the 
thermodynamic driving force for solute redistribution during the phase transformation and 
the diffusivity, respectively. If the latter vanishes faster than the former, the redistribution 
cannot be completely accomplished on the solid side of the interface, and trapping occurs. 
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We will now analyze this effect in more detail. 

Our starting point is Eq. (jU^ for the first-order diffusion field in the inner region. Without 
antitrapping current, its solution is 

f/i = f/i + ^ [1 + (1 - m] r rfe, (Al) 

which for the above choice of functions becomes 

f/i = f/i-t;„[l + (l-A;)f/o]r/. (A2) 

This solution, however, is not appropriate since it cannot be matched to the outer solution in 
the solid, which for a steady state is just a constant. The problem is that we have neglected 
terms in Eq. (jU^ that, for this solution, would not be small, which makes the calculation 
inconsistent. To see this, it is sufficient to remark that both the diffusion term (proportional 
to q) and the redistribution term (proportional to /i — 1) become exponentially small inside 
the solid. In contrast, for the above solution, the time derivative of Ui (equivalent to Vndr^ in 
the moving frame) gives a term of order e in the equation for f/, and hence becomes larger 
than the two mentioned terms far enough in the solid, which violates the counting of orders. 
In order to get a solution valid everywhere inside the solid, this term has to be included in 
the equation for Ui which becomes 

5,[g(0o)5,t/i] = ^ [1 + (1 - k)U,] 9,00 - ev^ - ^^0o) d,U, . (A3) 

By integrating once and using the boundary condition of vanishing current in the solid 
{q{4>o)d'qUi ^ for ?7 — > — oo), we find 

g(0o)5.f/i = y [1 + (1 - k)Uo] (00 - 1) - ev^ - ^^0o(o) d^\J^{i) di . (A4) 

For the sake of simplicity, let us first discuss the case = 1, in which the integral on the 
right-hand side is simply equal to V\{r]) — f/i(— oo). It can be seen immediately that this 
equation admits a solution that has the right limit, d^\J\ for ri —>■ —oo. We proceed 
by constructing an approximate solution by a matching procedure. First, remark that the 
left-hand side of Eq. ()A4|) is the product of two functions that vanish in the limit t] ^ —oo. 
Hence, it can be neglected in this limit, and the asymptotic solution is 

= U,i-oo) + l±i^=^(0o - 1) . (A5) 
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In contrast, in the region of the interface, the newly introduced term, being of order e, is 
small, which was precisely the reason to neglect it in the usual calculation. Therefore, in 
this region the solution of Eq. ()A2|) applies. Finally, a matching between the two solutions 
is found by searching the coordinate rj* where their slopes are equal, which, using the fact 
that 00 = — tanh(77/^/2), yields 

7]"" = -v^cosh"^ ^ . (A6) 



In the limit of small velocity f„, this simplifies to rj* = (1/^/2) ln(ef„/ ^/2). It can be checked 
that, in the matching region, the two terms (diffusion and time derivative) are of similar 
magnitude, which justifies the matching procedure. 

We have hence constructed an approximate solution, which is equal to the one obtained 
from the standard procedure for rj > rj*, and becomes a decaying exponential for rj < rj*. 
Evaluating the solvability integrals with this solution, we find, for example, 

F- = t/i(-oo) - f/i(0) = [1 + (1 - k)Uo] [l -2ln{evjV2)] . (A7) 

Similar terms appear also in the integral K. Using the identity eVn = {W/dQ){VndQ/D) = p, 
we find that the kinetic coefficient contains, in addition to the usual terms linear in p, 
corrections coming from eF" that scale as plnp. This constitutes, for small p, a logarithmic 
correction that makes convergence in p very slow. 

This calculation is an approximation, but the conclusion that there are nonlinear correc- 
tion terms is general, and can be easily interpreted: the anomalous kinetics occurs because 
solute can escape only from a region of size rj* behind the interface, and this size scales 
logarithmically with Vn (and hence p) in the limit of small p. In this limit, the case of 
arbitrary k can be easily treated and yields corrections of the form pln{kp). In the light of 
this conclusion, a physical sense can also be given to the condition used in the main body of 
the paper, namely that the "source function" must decay faster than the diffusivity: under 
this condition, all solute can escape the advancing interface. Note also that for a more 
realistic model in which the diffusivity becomes small but finite in the solid, the anomalous 
dependence of the kinetics on Vn stops for Vn < ?(+!) since then, again, all solute can escape. 
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APPENDIX B: DISCRETIZATION 



The phase-field and diffusion equations are discretized on a square grid of spacing Ax. 
We use standard finite-difference formulas, but a few details are worth mentioning. 
For the Laplacian of the phase field, we use the maximally isotropic discretization. 



which avoids the grid corrections to the anisotropy that are discussed in Ref. |13[. 

For the diffusion field U, we proceed by first calculating the current in each link, and then 
summing up all links around a site. On each link, the diffusion part, ju = —(pVU is calculated 
with the average of the phase-field according to ju = —{(pi+ij + 4>i,j){ui+ij — Uij)/{2Ax) 
for the X direction, and an analogous expression for the y direction. The most delicate part 
is the antitrapping current, jat = a{(j))W[l + (1 — k)U]hdt4>, where h = — V(/)/|V0| is the 
unit normal vector pointing into the liquid. We first evaluate the components of V0. The 
computation of the component parallel to the link is straightforward. As for the component 
perpendicular to a link, for a link along the x direction between sites and {i + 1, j) we 
use 



product a(0)[l + (1 — k)U]dt(p is then evaluated at the two end points of the link, and its 
average value multiplied with the appropriate component of n to obtain the current. 



[1] L.-Q. Chen, Annu. Rev. Mater. Res. 32, 113 (2002); W. J. Boettinger, J. A. Warren, C. Beck- 

ermann, and A. Karma, ibid. 32, 163 (2002). 
[2] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993). 

[3] W. Kurz and D. J. Fisher, Fundamentals of Solidification (Trans Tech, Aedermannsdorf, 
Switzerland, 1992). 

[4] J. Q. Broughton, G. H. Gilmer, and K. A. Jackson, Phys. Rev. Lett. 49, 1496 (1982). 
[5] L. V. Mikheev and A. A. Chernov, J. Cryst. Growth 112, 591 (1991). 



= -^{4>i+i,j + 4>i-i,j + 4>i,j+i + 4>i,j-i) 

+i(0i+ij+i + + 4>i+i,j+i + 0i-i,j-i) - 50, 



(Bl) 




49 



[6] J. J. Hoyt, M. Asta, and A. Karma, Inter. Science 10, 149 (2002). 

[7] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, Phys. Rev. A45, 7424 (1992); Phys. 

Rev. E47, 1893 (1993). 
[8] G. Caginalp and W. Xie, Phys. Rev. E48, 1897 (1993). 

[9] A. Karma, in Encyclopedia of Materials: Science and Technology, edited by K. H. J. Buschow, 
R. W. Cahn, M. C. Flemings, B. Ilschner, E. J. Kramer, S. Mahajan, Volume 7, Elsevier, 
Oxford, pp. 6873-86 (2001); in Thermodynamics, Microstructures and Plasticity, edited by A. 
Finel et al., Kluwer Academic Pubhshers, North Holland, pp. 65-89 (2003). 

[10] J. Tiaden, B. Nestlcr, H.-J. Diepcrs, and 1. Stcinbach, Physica D 115, 73 (1998). 

[11] S. G. Kim, W. T. Kim, and T. Suzuki, Phys. Rev. E 60, 7186 (1999). 

[12] J. Bragard, A. Karma, Y. H. Lee, M. Plapp, Interface Science 10, 121 (2002). 

[13] A. Karma and W.-J. Rappel, Phys. Rev. E 53, R3017 (1996); Phys. Rev. E 57, 4323 (1998). 

[14] N. Provatas, N. Goldenfeld, and J. Dantzig, Phys. Rev. Lett. 80, 3308 (1998); J. Comp. Phys. 
148, 265 (1999). 

[15] M. Plapp and A. Karma, Phys. Rev. Lett. 84, 1740 (2000); J. Comp. Phys. 165, 592 (2000). 

[16] A. Karma, Y. H. Lee, and M. Plapp, Phys. Rev. E61 3996 (1999). 

[17] J. A. Warren and W. J. Boettinger, Acta Metall. Mater. 43, 689 (1995). 

[18] R. F. Almgrcn, SIAM J. Appl. Math 59, 2086 (1999). 
[19] A. Karma, Phys. Rev. Lett. 87, 115701 (2001). 

[20] W. Kurz and D. J. Fisher, Fundamentals of solidification (Trans Tech, Aedermannsdorf, 
Switzerland, 1992). 

[21] M. J. Aziz, J. Appl. Phys. 53, 1158 (1982); N. A. Ahmad, A. A. Wheeler, W. J. Boettinger, 

and G. B. McFadden, Phys. Rev. E 58, 3436 (1998). 
[22] K. R. Elder, M. Grant, N. Provatas, and J. M. Kosterlitz, Phys. Rev. E 64, 021604 (2001). 
[23] G. B. McFadden, A. A. Wheeler, and D. M. Anderson, Physica D 144, 154 (2000). 
[24] W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 35, 444 (1964). 

[25] S-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Corieh, R.J. Braun, and G.B. 

McFadden, Physica D 69, 189 (1993). 
[26] R. Folch, J. Casademunt, A. Hernandez-Machado, and L. Rami'rez-Piscina, Phys. Rev. E 60, 

1724 (1999); 60, 1734 (1999). 
[27] M. Georgelin and A. Pocheau, Phys. Rev. E 57, 3189 (1998); Eur. Phys. J B 4, 169 (1998). 

50 



[28] R. Folch and M. Plapp, Phys. Rev. E 68, 010602R (2003). 

[29] J. C. Ramirez, C. Beckermann, A. Karma, and H.-J. Diapers, Phys. Rev. E 69, 051607 (2004). 

[30] Y.-T. Kim, N. Goldenfeld, and J. Dantzig, Phys. Rev. E 62, 2471 (2000). 

[31] G. Amberg, Phys. Rev. Lett. 91, 265505 (2003). 



51 



